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ABSTRACT 


A comparison of numerical methods utilized by the finite 
element technique for solving a nonlinear nuclear reactor 
dynamics problem was conducted. Using the Crank-Nicolson, 
DVOGER (Gear) and Implicit Gear methods, the results showed 
tne implicit to be the superior method investigated. This 
is based on the fact that all three methods yielded the same 
steady state solutions; but, the Implicit Gear method used 
significantly less CPU time and comparable storage to Crank- 
Nicolson. This was particularly apparent as the degrees of 
freedom were increased. In addition, the transient solution 
in all cases was better than that obtained in Crank-Nicolson 
and compared favorably to that of Gear's method. 

The other noteworthy result was in the effect of the 
error criterion on solution. It was shown that for a range 
of error from 1 i to 1.0, the steady state solution value 
remained the same. This results in a significant reduction 
in computer processing time since the time required decreases 


substantially as the error conditions imposed are relaxed. 
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I. INTRODUCTION 


A. PURPOSE 

This research project has been undertaken to compare 
several numerical methods of solving a nonlinear nuclear 
reactor dynamics problem. Three methods have been investi- 
gated in this thesis, including the Crank-Nicolson, the 
DVOGER (Gear) and the Implicit Gear methods of solution. 

A nuclear reactor dynamics problem with temperature- 
dependent feedback, when it entails either a non-homogeneous 
or multi-region reactor, results in a nonlinear field equa- 
tion in space and time. This problem does not lend itself 
to solution by analytical means [5, 6]. However, when the 
physical and neutronic properties of the problem are known, 
a model can be formulated using the finite element method 
which will yield the transient and steady state flux solu- 
tions. In particular, the partial differential equations 


investigated were of the form 
at = bY? Ws cy aw" (1) 


meuietc 42,0, ©, dud W are constants and W (rv, z, t) 1s the 
flux. The finite element method reduces Equation (1) to 


the system of ordinary differential equations 


N N 


sy Of the. Go fy agcepy. — rot 2 
Deere 2 seer y. feyeneee SUNUT PETE 2°¢2) 


when the nonlinear term is linearized. Nguyen and Salinas [5] 
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and more recently Olsen [6] discuss the problem and methods 
of solution. 

In this work, a comparison of three numerical integration 
schemes for Equation (2) was made. The comparisons will be 
made on several bases: computer storage requirements, com- 
puter processing time, rapidity with which solution was 
obtained and the relative accuracy of the solution. 

At steady state, it is expected that all three methods 
will provide the same solution value. This is due to the 
fact that at steady state the time derivative of the flux 
(wv 7; Bauation (2), is Zero. 

In order to explore the relative value of the various 
methods, the model has been discretized into finite element 
grids of various degrees of freedom (DOF). The effect on 
the solution by finer discretizations of the finite element 
has been investigated to determine if the various methods of 
wetueion dre simtlariy mituenced ‘to provide a’ ‘better ‘solution 
fer a finer mesh. To have a method of solution relatively 
independent of mesh size would greatly reduce computer storage 
and processing time requirements since a larger grid with 
fewer elements could be utilized. 

In order to test the flexibility of the various equation 
solvers, an initial disturbance was introduced at different 
points in the model. This provided both a check of the 
ability of the method to accept a random disturbance as well 
as information on how rapidly a solution is obtained with a 


particular disturbance input. Two nodal points of the system 
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were considered, a point at the origin and a point on the 
core-reflector interface. This was done to determine if the 
tracking ability was consistent throughout the model. 

Additional comparisons investigated include the effects 
of the convergence criterion on the solution for all methods 
and the effects of the size of thie time step on the Crank- 
Nicolson method. 

Modifications have been made to the programs provided in 
Ref. [6] which was the basis of this research project. In 
that work, core properties were arbitrarily assigned »to the 
reflector elements at the interface. This occurred because 
the material and nuclear properties were provided at the nodal 
paanes tather than. by elements. .The properties were mtro- 
duced in this manner as a means of computer storage reduction. 
Regardless of the mesh size, there were always fewer nodes 
than elements. In this project, all properties were provided 
on an element basis to eliminate this discrepancy while at 
the same time sacrificing the minor storage savings it repre- 


sented. 
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II. DATA GENERATORS 


A. GENERAL 

The data generators formulated in this work are useable 
only for regular rectangular grids (Figures 1-3). Having 
selected the number of horizontal and vertical points for 
the discretization of the model, the data generators will 
provide the numbering of each nodal point, the horizontal 
and vertical position of each node and the nodal neighbors 
of each node. 

In addition, the outer boundary nodes, at which the 
dynamic flux is zero, will be numbered last in order to 
reduce the number of equations required by the particular 
method of solution being used. This will substantially 
@eeaecase the Storage requirements..." For ‘example, in a one 
hundred thirty-two node discretization, only one hundred ten 


equations will be solved. 


B. PROPERTY INPUTS 
1. Purpose 

A simple data generator has been provided which will 
produce a data deck containing the physical properties of 
the reactor core and reflector in the format required by the 
Crank-Nicolson and DVOGER (Gear) methods. These properties 
are neutron velocity, neutron diffusion length, neutron 
absorption cross-section, neutron fission cross-section and 


reactivity temperature coefficient. 
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This generator will accommodate any size mesh and 
will process an indefinite number of grids simultaneously. 

The user must provide this routine with a data deck 
which contains the total number of elements in the grid and 
the type of each element, eore or reflector. This has been 
simplified oy the designation of all core elements by ITYPE=0 
and reflector elements by ITYPE=1. 

2. Programming Notes 

a. The property values in the two regions must be 
provided as an integral part of the program. 

b. The type of element, core or reflector, must be 
provided. 

c. The routine will process an unlimited number of 
models. It, therefore, must be halted when all data has been 
utilized. This is done by specifying the final value of the 
number of elements in an IF-STOP statement. 

as. “Parameters 

a. NEL - number of elements in the discretized model 

betters = type of element; ITYPE=0 is a core element 
and ITYPE=1 is a reflector element. 

c. D - vector containing the diffusion length of the 
core and reflector elements. 

d. SGA - vector containing the absorption cross- 
section of the core and reflector elements. 

e. SGF - vector containing the fission cross-section 
of the core and reflector elements. 


i) = Vector containing the neutron, velocity. 
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g. ALPHA - vector containing the reactivity 
temperature coefficient. 


C. NODAL POINT COORDINATES AND ELEMENT 
NODAL POINT CONNECTIVITY 


i. Purpose 

| Nodal point positioning in the model is readily 
obtained in data deck form from this routine. Once the model 
to be investigated has been discretized into the finite ele- 
ment grid desired, the user provides the number of vertical 
and horizontal nodal points and their dimensional position 
along the axes. By the use of a nested loop, the nodal points 
are consecutively numbered and assigned the proper coordinate 
dimensions. 

The element connectivity, that is, the nodal points 
for each element, is also resolved by the use of a nested 
loop and several counters. Each iteration will yield two 
elements with their respective nodal point boundaries. This 
will continue until the nodal points for each element have 
been computed. 

The output of this routine will consist of two data 
decks. The first will provide the nodal point with its ver- 
tical and horizontal coordinates. The second will yield the 
element nodal point connectivity. 

Omitted from the element connectivity data deck, but 
required by the thesis program, is the type of element, core 


or reflector. This must be added to the deck individually. 
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Example of output (See Figure 1): 


nodal point coordinates 
nodal point horizontal position vertical position 


i 0.0 0.0 


element connectivity 
element nodal point connectivity type element* 


1 i 9 10 0 


*Note: must be added manually to data card. 


Fe Programming Notes 


a. The maximum number of vertical and horizontal 
nodal points in the discretized model must be provided. 

Be The herazontal and vertical positions of the 
discretization must be provided. | 

c. The boundary points of the model will be numbered 
such that these nodal points are the last in the numerical 
sequence. 

du» (fhe. routine.isS written such,that,it,will.process 
an indefinite number of models. Therefore, it must be halted 
when all the data has been utilized. This will be accomplished 
by specifying the maximum number of vertical points in the 
model in an IF-STOP statement. 

5, (Parameters 

a. NH - maximum number of horizontal nodal points 
in the discretized model. 

b. NV - maximum number of vertical nodal points in 


the discretized model. 
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c. X - vector containing the horizontal positions 
of the NH points. 

a- Y= vector containing the vertical positions of 
the NV points. 

e. SYSNOD - vector providing the total number of 
nodal points. 

£-.. Rk > Vector providing the radial. position of the 
nodal points. 

g- .4 - wector providing the. vertical position of 
the nodal points. 

h. ELNOD - array providing the nodal point boundaries 
for each element. 

i. NEL - the total number of elements in the discre- 


tized model. 


D. NODAL NEIGHBOR CONNECTIVITY 
1. Purpose 

This’ routine produces a data deck for the Crank- 
Nicolson and Implicit Gear methods containing each nodal 
point and the nodal points connected to it by the finite 
element discretization. 

The regular rectangular grids are such that no more 
than six nodal points contribute to any one; therefore, an 
Nx7 array can be formulated for nodal neighbor connectivity. 
For example, nodal point 26 in the 112 element grid would be 
stored as follows (Figure 4): 


26 Ney 18 ah 5D 34 25 
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As written, zeros are inputted whenever a nodal point does 
not have six neighbors to maintain the symmetry of the matrix. 
This could, however, be compacted further by reading the data 
as a single vector with a starting point counter to indicate 
when the next nodal point has been reached in the sequence. 
This would eliminate the need for the zeros to be supplied. 

This routine will provide a data deck for any 
rectangular grid and for an indefinite number of grids. 

2. Algorithm 

This algorithm provides a dense, compacted matrix 
which reduces the storage size required by the main program. 
The use of a finite element method allows a certain amount of 
Eempacting, i.¢., the banded matrix. This is a function of 
the finite element size, that is, the eC eaeee the number of 
nodal points, the larger the band. The size of the band is 
determined by the largest difference between nodal neighbors 
plus one. For example, in Figure 1, node ll has a maximum 
difference of 20 - 2 = 18, while at node 16 the difference 
= 45 — 7 = 50 Which 15 the Largest difference in the forty- 
five node system. As the number of nodal points increases, 
this maximum band width also increases. In Figure 2, the 
Maximum difference is at node 16 (67 - 7 = 60}; and, in 
Figure 3, the maximum difference of 124 - 10 = 114 exists 
Beemoae 22. Thereétore, although the size of the system is 


reduced from NxN to Nxq, the size is not constant and may 
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not be small. In general, if one used banded storage, the 
numbering would proceed consecutively in the direction of 
fewest nodes. However, this would require the identification 
of those nodes on the boundary since they are of constant 
value and do not require integration. 

In the particular scheme of discretization utilized 
in this project, no nodal point has more than six nodal point 
contributors. This allows the matrix to be compacted further 
from Nxq to Nx7. 

3. Programming Notes 

a. The total number of nodal points, the maximum 
number of nodal point contributors and the maximum number of 
horizontal and vertical nodal points must be provided. 

b.. The routine will satisfy any rectangular grid. 
However, it must be instructed when a sufficient number of 
nodal points have been generated. This is accomplished by 
using the total number of nodal points desired in an IF-GO TO 
statement. 

c. This routine requires a positive means of stopping. 
This is accomplished by the use of the maximum number of 
vertical nodal points in the last data set in an IF-STOP 
statement. 

dz) ime output of this routine.will place the central 
nodal point first with the contributors following. 

4. Parameters 
a. NV - number of vertical nodal points in the model. 


b. NH - number of horizontal nodal points in the model. 
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c. NVH - total number of nodal points in the model. 
NVH = NV x NH. 


d. LCON - maximum number of contributing nodal points. 


e. MNOD - an array, NVH x LCON, providing the central 


and contributing nodal points. 
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III. CRANK-NICOLSON METHOD OF SOLUTION 


A. DESCRIPTION 

The Crank-Nicolson formulation is a numerical method 
SFiginally presented by J.*Crank and P.. Nicolson in 1947. 
Crank-Nicolson, being an implicit method, does not require 
the inverse of the matrix to be calculated. Therefore, 
advantage is taken of the sparse matrix inherently provided 
by the finite element method. If the inverse of the sparse 
mectex 35) formed;) ve) wili: be. full. 

The Crank-Nicolson method will solve the set of linear 
differential equations represented by 
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Crank-Nicolson, after some algebra, yields 


pee a | ee ey (4) 
At 2 

The implicit system (4) may be solved by an iterative process; 
in this case, the Gauss-Seidel method is used to solve the 
set of simultaneous algebraic equations (4) formulated. Of 
all the stable numerical methods in the case of single step 
implicit, Crank-Nicolson has the smallest truncation error [2]. 
The Gauss-Seidel method of iteration is continuously using the 
newest solution values allowing convergence to the solution 


to be the most rapid. 


Ze 


? ¥ 
7x) 


1 ; rr > 
DE TUIOe YO conrel tent 0 AMAR iT 
| worremand, 


soiseilomro3 nozloztheanez9 if 


inst) .L yd betaeeorg visa 


tem tiaiiqnd ne gated soetoath~aes 
isluoige sd as sees st eiy- to eerevink i 


‘€ = . : ae 


irtem Se7raqe #d2 Yo pete? 22 ogsthsv 
oe 
oft ti  hoittem snomele otini? saz a 


Iliw 3% .bDsmr6i ef a1t38 


a 


5) ‘f2 f . P + a Hie fA ri: need sit 


v . + 
Sisters wissupes Lalineyei 


4 sane 45 . foelosiv-danar 


. 

’ ls, <a ~ 39 

— —— | a ee en “ar heed 
| 


itepo otetdegis gucensilumte te a 
af siuotyem feo rremun sldat2 nit ta 


ent sor phates bit: iaB%s0 abgeh iqime 


rghrerea | 10. borden reblepeeue 


a. 


. 63 erevans gtlwoble -egulev ents = 


bik r08 
" 


- Pe, 


B. .EFFECT OF DEGREES OF FREEDOM 

As might be expected, as the grid size became finer, the 
computer core requirements increase. In addition, more time 
is required for the problem to reach a steady state solution. 
At the same time, as shown in Figures 5-7, the solution values 
converged, as the mesh size became finer, to a more accurate 
solution. 

As shown in Table I and Figure 8, storage and processing 
time, that is, CPU time to reach a steady state solution, in- 
creased markedly between a forty-five node grid and a one 
hundred thirty-two node grid. This yielded a fifty-nine 
percent increase in storage and better than a five hundred 
percent increase in processing time. At the same time there 
is an improvement in the solution of only four and one-half 
percent. 

Comparison to a seventy-two node grid yielded far more 
satisfying results. At seventy-two nodes, there is an eighty 
percent increase in CPU time and a sixteen percent increase 
in storage requirements while obtaining a two and one-half 
percent improvement in the solution. 

Similar results were obtained in problem time to solution. 
For one hundred thirty-two nodes, the time more than doubles; 
while for seventy-two nodes, the time was increased by twenty 


percent. 


C. EFFECT OF ERROR CRITERION 
The error criterion used throughout this thesis informs 


the integration routine when it has achieved sufficient 
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agreement between iterates (i.e., when it can stop and go 
to the next iteration step). 

In Crank-Nicolson, the error:criterion simply consisted 
of the difference in successive Gauss-Seidel iterates divided 
by the current iterative value. 

The error criterion was varied from jos to Vie to 
observe the effects on the solution and computer requirements. 
For Crank-Nicolson the steady state solution value remained 
the same regardless of the error criterion. This demonstrated 
that the extra time required to satisfy a tighter error cri- 
terion at each eee was not necessary to reach a satis- 
factory steady state solution. As was expected, the computer 
processing time did increase significantly as thevercon 
criterion was made smaller. An unusual result occurred when 
4 


the error criterion was set to 10 This particular value 


resulted in a processing time less than. that required for 
a, It would appear that this might have been the result 
of two contributing factors. With the closeness of the error, 
each time an iteration was performed, it was using a better 
solution, and each new time step also had a better solution 
value. Although there were many more time steps required 

for this error criterion, the steady state value was rapidly 
approached because fewer iterations were required at each new 
time. These results are shown in Table II. These results 


= ‘3 - -4 
for an error criterion between 10 1 and 10 


do not imply 
that this would be the result obtained for all problems of 


this type. 
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Be EFFECT OF INCREASE IN TIME STEP SIZE 

In this method, the time step can be increased or not 
depending on the solution. If a solution is not obtained 
with a particular time step, that same time step is reduced; 
and the routine attempts to attain solution with the smaller 
time step. This continues until either a satisfactory solu- 
tion value is obtained or the built in default value is 
reached which terminates the routine. When a solution is 
attained, the routine compares the number of iterations 
required to reach that value to that of the previous time 
step and to a specified number of iterations. If the current 
number of iterations is less than either of those numbers, 
the program increases the time step by an input constant 
value. 

The change made to the initial time step as the steady 
state was approached was the critical value if the solution 
was to converge to a steady state. It appeared that Crank- 
Nicolson was particularly sensitive to the amount the time 
step was increased as the solution was approached (Figure 9). 
The method was unable to recover if, when nearing the knee of 
the curve, the size of the increase was such that the steady 
state value was exceeded. Once the solution was passed, the 
results indicated a diverging oscillation about the steady 
Stace value. 

Several values of the initial time step were utilized in 
all grids to attempt to locate the steady state. The results 


Of these trials are provided in Table III. It is noteworthy 
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that if the increase was greater than 1.2 times the initial 
step size, the steady state would never be achieved. Whether 
Emes would be the case for all grids is not obvious. This is 
particularly evident in view of the results obtained by 

Olsen [6], who utilized an increase of 1.5 times the initial 
step for a thirty-eight node grid. It would, therefore, seem 
that a trial and error approach would be necessary until an 


increase in step size resulted in a steady state solution. 
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IV. DVOGER (GEAR) METHOD OF SOLUTION 


A. DESCRIPTION 

The DVOGER (Gear) routine is an IMSL library routine 
which integrates a system of explicit first order differ- 
ential equations. Within this library routine is the capa- 
bility of solving both stiff and nonstiff systems by selection 
of an indicator. In solving the nonstiff system, the Adams 
predictor-corrector is utilized. For the stiff system, Gear's 
predictor-corrector is used. Gear's method computes the 
Jacobian for the system of ordinary differential equations 


mm order £0 optimize the time step for each integration. 


B. EFFECT OF DEGREES OF FREEDOM 

Gear's method placed a much greater demand on the computer 
in terms of processing time and storage than the other methods 
considered. This was primarily due to three causes: 1) the 
calculation of the Jacobian, 2) the transformation of Equation 
om <e explicit Lom, did 3) the initial ‘absolute value of 
the solution. The first two items require a substantial 
increase in both CPU time and computer storage. Several 
variances were implemented in this program to determine if 
better use could be made of the method to make it competitive 
with Crank-Nicolson or the Implicit Gear methods. 

When the system was considered nonstiff, the storage 
requirements decreased since the Jacobian was not calculated; 


but, the processing time increased by one order of magnitude 


Zi 


it a 
ys PP eS 
~ ‘J rt] H 2 
~ =< NL! F tt - 
rad a ‘ 
; oo a, 
Ibe | ORTHM ¢%. AID 3D) Axo0vG  . Vr” 
ne —— - | tttndemetineemee ae a eo 
? rc ' uu .> . 
c YOTT Dea 
a4 ) ) Rav 
2o7eTgOIal Ankaw 
a | 
- Y Pad7im | .2nolssuos: [eisne 
> A 
ad gntvloe to \itiie 
. Me 
t stoas 
a 5S FTto-: rot > ibe 1g 
ji xt ‘1 50D>-t0s> Ebsse 
£5 bi 4 fy 70? nsidosat 
woh TO 
5 ; 3 7 70 
7 & ¢ 19] 


fy (f Jpeiteset Lo moisaiugigss 


via sais UID diod ni sesetont 


> ha 


‘teqonmeldal etew seonnigs 


. . ‘ haa A Bilins | reszed: 


: 2254 ir 30) weelosty- Anexd dake 


when compared to the stiff system treatment due to the small 
time steps taken to solution. The routine calls for an 
initial absolute solution value of one to be supplied. This, 
however, adversely affects the progress of the problem since 
this value of one is initially compared to values of right 
fats Jimits the time step taken by the routine initially and 
continues to affect the progress until the steady state solu- 


tion is neared. When an initial value of 1024 


was utilized, 
the processing time was decreased by forty percent. This was 
due to the fact that larger time steps were allowed from the 
outset, since the previous solution value (in this case the 
initial value) was comparable to the value obtained in the 
first calculation. 

At best, when utilizing Gear's method, the processing 
times were two orders of magnitude greater than that required 
for the Implicit Gear method and twenty times the Crank- 
Nicolson processing times for the 132 DOF system, as shown in 
Table I. Core requirements were also increased significantly, 
particularly at one hundred thirty-two degrees of freedom. 
For this grid, DVOGER (Gear) was approximately three times 
Crank-Nicolson and the Implicit Gear core requirements. 

Gear's method provided the same peeaay, eaee solution 
values as the other two methods, and the transient curves 
were very similar to those obtained in the Implicit Gear 


method (Figures 10-15). 
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C. EFFECT OF ERROR CRITERION 

In Gear's method, the time step size is adjusted so that 
the single step error estimate divided by the previous maximum 
solution value is less than the error criterion in the 
Euclidean norm. The single step error estimate is a multiple 
of the difference between the predicted and corrected values 
of the variable. 

fae -+error criterion, fable II; was.varied,for.this,method 
as it was for the others. The results were the same; as the 
criterion was ee CPU time increased. In this case, 
however, the time became excessive very rapidly, more than an 


hour for 10> and almost four hours for 10 * for the model 
with 45 DOF. Similarly, the solution values were the same 


regardless of the criterion utilized. 
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V. IMPLICIT GEAR METHOD OF SOLUTION 


A. DESCRIPTION 

This numerical method is particularly useful for the 
solution of large, sparse systems of implicit, stiff differ- 
ential equations [4]. In contrast to Gear's method, the 
Implicit Gear method eliminates the need of an exact NxN 
Jacobian, but rather requires only an exact Nx7 Jacobian. 
This is due to the fact that Equations (3) are handled directly, 
thereby retaining the sparseness of the matrix. 

Gear's predictor-corrector method and the compact matrix 
makes use of storage because the implicit Equations (3) are 
handled directly. This results in very efficient use of the 
computer both in terms of storage and processing time require- 
ments. 

The user ts required to provide a subroutine that evalu- 
ates the system of equations being investigated, as well as, 
for efficiency, a subroutine to evaluate the Jacobian of the 
system of equations. 

This program is a modification by Franke [4] to the 
routine DFASUB [7]. One of the major changes to the routine 
is in the treatment of the error. In the Implicit Gear 
method, instead of using the Euclidean norm, the root mean 
Square norm is utilized. This is no more than the Euclidean 
norm divided by the square root of the number of components. 
In addition, the maximum value of the component is updated 


before the norm of the relative error is computed. 
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B. EFFECT OF DEGREES OF FREEDOM 

The steady state values obtained in this method were the 
same as those obtained in Crank-Nicolson and DVOGER (Gear) 
as shown in Figures 11-18. This result does not imply that 
these methods will always give identical steady state solu- 
tions. Certainly they do give different transient solutions. 
The major advantages of this method are 1) the computer 
processing time was reduced as the grid size became finer 
and 2) the storage requirements are slightly greater than 
that required in the Crank-Nicolson method due to the compu- 
tation of the Jacobian and the storage of up to seven past 
Seiutions which control the size of the time step... Implicit 
Gear requires about Nx25 more storage locations than Crank- 
Nicolson. This amounts to about only thirteen thousand bytes 
for the one hundred thirty-two degrees of freedom system 
(single precision). The magnitudes of the increase in pro- 
cessing time dropped dramatically in this method. For the 
same initial disturbance, the time increased eighty percent 
for the one hundred thirty-two node grid and twenty-four 
percent for the seventy-two degrees of freedom. In addition, 
Implicit Gear gave a more accurate transient solution. The 
results of this method are shown in Table I and comparison 


to the other methods will be made in Chapter VI. 


> EPPECT OF ERROR CRITERION 
im imptiert Gear, the error criterion is defined as in 


DVOGER (Gear) except that the root mean square of the Euclidean 
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norm and the updated maximum solution value are used in the 
computation. 

The error criterion was varied substantially to determine 
what effect it had on computer storage and time requirements, 
maeewell as ats effect on the accurdcy of solution. As shown 
in Table II and Figures 19-20, a wide range of convergence, 
1.0 to to", was investigated with interesting results. 
Expectedly, the processing time did increase as a tighter 
error criterion was required. However, the criterion had 
little effect on the track through the transient solution for 
values of less than ioe, When using an error of one-half 
and one, the transient solution varied substantially. The 
same steady state solution was obtained, although at a later 
point in problem time as shown in Figure 19. For the stiirer 
error criterion requirements, the processing time necessary 
is doubled by utilizing a criterion of to * instead of a value 

1 


greater than 10 ~. This, certainly, did not appear to warrant 


the closer tolerance level. 
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A. COMPARISON OF TIME AND STORAGE REQUIREMENTS 

With all three methods yielding the same steady state 
solution (Figures 10-18), the comparison of methods became 
one of time and storage requirements placed on the computer 
rather than one based on which provided the best solution. 
In the transient stage, Gear tracked slightly better than 
Implicit Gear, but both Gear and Implicit Gear tracked better 
than Crank-Nicolson (See Figures 10-18). For the three systems 
of equations utilized, the Implicit Gear method performed 
Significantly better than either Crank-Nicolson or DVOGER 
(Gear) in CPU time and was comparable to Crank-Nicolson and 
superior to Gear in storage requirements (Table I). This 
becomes more significant as the number of degrees of freedom 
increase. The Implicit Gear is less sensitive to the increase 
in size of the system of equations as shown in Figure 8. 

Comparing the forty-five and the one hundred thirty-two 
degrees of freedom, the core requirements increased slightly 
and the time doubled for the Implicit Gear method. In con- 
trast, for Crank-Nicolson, the time increased by six times 
aid the’ core by LOOK bytes; and, for DVOGER (Gear), there 
was a ten times and 500K byte increase in time and storage 


respectively. 
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B. DEGREE OF FREEDOM EFFECT ON SOLUTION 

Varying the degrees of freedom resulted in a better 
solution value as shown in Figures 5-7 and 21-26. There was 
a four and one-half percent better solution obtained for the 
one hundred thirty-two degrees of freedom compared to the 
forty-five degrees of freedom. This better solution is very 
costly both in terms of computer core requirements and CPU 


time (Figure 8). 


C. ERROR CRITERION 

Error criterion was varied for all methods to determine 
the effect on solution and computer. As might be expected, 
the computer processing time increased for all methods as 
more rigid tolerance was imposed (Table II). However, although 
the criterion was made very close, there was no appreciable 
effect on the transient solution until the error criterion 
was greater than io and no effect on the steady state 
Seigeton. this is significant in that, as stated above, the 
processing time increases greatly as the criterion is tight- 
ened. This shows that for the particular problem considered 
here, some close error criteria yield no advantage, only the 


disadvantage of requiring more time in the computer. 


D. INITIAL DISTURBANCES 

The input of various initial disturbances (central, skewed 
at the core-reflector interface and uniform throughout the 
core) yielded the same steady state solution value. The 


transient solution varied substantially due to the input 
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position of the disturbance. In most cases, due to the 
Magnitude and scope of the initial disturbance, the uniform 
disturbance required the least processing time to reach the 


steady state value. 


E. COMPARISON OF SOLUTION TRACKING 

Two nodal points were investigated graphically for each 
disturbance, integration method and DOF, as shown in Figures 
a7, 10-18 and 21-26. At these test points, the methods 
performed similarly in all cases in both the transient and 


steady state solutions. 
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Vit: TEST PROBLEMS 


A. PHYSICAL MODEL 

A cylindrical reactor with the dimensions and properties 
Given in Table IV was used as a model. A radial slice of 
this model was discretized by various finite element grids 


as shown in Figures 1 =3. 


B. COMPUTER PROCESSING CONSIDERATIONS 

The programs presented in this thesis, Appendix A, were 
written in the FORTRAN IV language and all computer runs, 
Table V, processed on the IBM 360/67 computer using the 
FORTRAN 'G* compiler. It is expected that the FORTRAN 'H' 
compiler would have supplied results similar to those obtained 
with the FORTRAN 'G' compiler. This was not done because the 
FORTRAN 'H' compiler requires 350K bytes as a minimum core 
requirement. The majority of the runs performed in this 
thesis required significantly less than 300K storage. Single 
precision (six to seven significant digits) was used through- 
out this research. As has been shown [6], double precision 
solutions are at variance by less than 0.01 percent with 


single precision solutions. 


C. PROBLEM ANALYSIS 
The techniques used by Olsen [6], modified as described, 
and the Implicit Gear method [4] were utilized to solve the 


problem delineated in Ref. [5]. The problem was approached 
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using finite element discretizations with forty-five, seventy- 
two and one hundred thirty-two degrees of freedom and providing 
an initial disturbance. Three initial disturbances were pro- 
vided as follows: central at the core center (R = 0 cm., 

Z = 0 cm.); uniform throughout the core; and, a skewed dis- 
meurbance at the core-reflector interface (R = 60 cm., Z = 0 


em.). 


D. PROGRAM USAGE 

In order to properly utilize the routines presented in 
this thesis, several points must be considered. 

In the Crank-Nicolson routine, there was difficulty in 
obtaining a steady state solution. The ultimate cause was 
the size with which the previous time step increases. This 
requires a trial and error approach for the particular grid 
Size. For this project, an increase of ten or twenty percent 
yields satisfactory results; but, anything larger results in 
a divergent oscillation about the steady state value. 

In all three methods, the dimensioning should be set by 
the particular problem being solved, i.e., determined by the 
degrees of freedom. If this is done, most efficient use will 
be made of the computer, and the user will experience better 
turn around time on the equipment. 

In the use of the Implicit Gear method, the data cards 
for nodal neighbor connectivity must have the contributing 
nodes listed consecutively, with the zeros, used for nodes 
not having six neighbors, being last on the cards. For example, 


in the one hundred thirty-two degrees of freedom system 
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(Figure 3), the nodal neighbor connectivity for node 122 


would be written as 


£22 2S 41 0 0 0 OF 


The data deck arrangement for all three methods is given in 


manieVi.< 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

It has been shown in this research project that of the 
three methods investigated in the solution of the test problem 
that the Implicit Gear performed in a superior manner in CPU 
time requirements and was comparable to the storage require- 
ments of Crank-Nicolson. In addition, Implicit Gear was the 
least sensitive to the change in the degrees of freedom. 

The method used had no effect on the steady state solution 
value obtained as all three yielded the same results for the 
same DOF. In the transient solution, Implicit Gear conformed 
very closely to DVOGER (Gear). As shown in Figures 10-15, as 
the DOF increased, the transient solutions moved in the direc- 
tion of the transient solution of Gear's method. 

As the DOF in the discretized finite element was increased, 
a better solution was obtained. This, however, was counter- 
getea by the fact that for the small increase in accuracy 
obtained by the finer mesh, the time and storage need of the 
computer increased significantly. 

A very important Fesult of this project was the BE ReCE of 
varying the error requirements of the problem. This yielded 
essentially no change in the accuracy of the transient solu- 
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to 10 ° and no variance in the steady 
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methods of integration, the error criterion selected would 
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be a function of the solution desired. If the user is only 
concerned with the steady state solution, the tolerance could 
be relaxed to 1.0. However, if the transient solution is 
needed, an error of ie appears to be the least rigid value 


which will still provide a satisfactory track to steady state. 


B. RECOMMENDATIONS 

It would be of value to investigate further the Implicit 
Gear method through the use of other problems solving a non- 
linear system of ordinary differential equations. In addition, 
larger degrees of freedom should be attempted to evaluate any 


limitations that may exist in this method of integration. 
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TABLE I 


COMPARISON OF SOLUTION METHODS 


METHOD CORE (K) PROCESSING TIME* (min.) 
Central Skewed Uniform 
45 Nodes 
CRANKO 158 67 7S 3357 
DVOGER (GEAR) 148 L728 19..5 255 
IMPLICIT GEAR 124 1.05 1.00 1320 
72 Nodes 
CRANKO 184 5,07 5.49 2.58 
DVOGER (GEAR) 244 89.5 114. d 59.2.5 
IMPLICIT GEAR 124 ies 12350 13.25 
132 Nodes 
CRANKO 252 ees oie G4 
DVOGER (GEAR) 610 >400 >400 >400 
IMPLICIT GEAR 124 1.9 20 2.0 


Processing time is that required tod reach a steady state 
solution. 
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TABLE II 


EFFECTS OF SOLUTION ERROR CRITERION 


45 Nodes 
SOLUTION ERROR PROCESSING TIME TIME AT SOLUTION 
CRITERION (min. ) (sec.) 


CRANK-NICOLSON 


0.1 1.8 5569 2 40 
0.01 5.5 S70 x 10” 
0.001 14.2 5°57 x 10° 
0.0001 0.83 5.87 x 0 
IMPLICIT GEAR 
1.0 0.99 7% 56 4 LOre 
0.5 0.92 Toe ee alee 
0.1 1.05 Tee x aVOlue 
0.01 Tao fade Mae 
0.001 1.45 Ga56 x AO ne 
0.0001 1.58 ee44 x for 
DVOGER (GEAR) 
0.1 18.6 Bet eo 
0.01 255 579 x 10° 
0.001 78.8 aos 2 0 
0.0001 >300 = 
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1000! 
1000.0. 


onthe 
2.000 
‘an 
“10.0. © 
100,04 
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TABLE III 


Bepecits OF VARYING THE INCREASE OF THE INITIAL TIME STEP 
IN CRANK-NICOLSON METHOD 


INCREASE CORE PROCESSING TIME SOLUTION TIME 


(min. ) (sec.) 
45 NODES 
1.1 158K 4.99 3.18 x 10°° 
1.2 158K es eo te 
is 184K 7.42 OSCILLATING 
72 NODES 
1.1 252K 5S Ase Oe” 
12 184K 3.05 Bs x 10 
Oke 252K 205 OSCILLATING 
1.4 252K 9.1 OSCILLATING 
1.5 184K 621 OSCILLATING 
132 NODES 
bt 252K S54 5.58 x 10° 
1.2 252K 8.4 8.03 x 10° 
iJ5 184K 2.6 OSCILLATING 
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SYMBOL 


= fF WD Ww 


ac 


ar 


COMPUTER 
SYMBOL 


c-N, UP 7 IG 


R 

R 

Z 

Z 
V/VELOCT 
D/DSUBF 


D/DSUBC 


SGA/SGMAAF 


SGA/SGMAAC 


ZNU 


SGF/SGMAF 


SGF/ - 


FISFAC 
HBAR 


ALPHA 


TABLE IV 


PHYSICAL CONSTANTS 


DEFINITION 


Total radius 

Core radius 

Core height 

Total height. - 
Neutron velocity 
Neutron diffusion 


coefficient (core) 


Neutron diffusion 
coefficient 
(reflector) 


Neutron absorption 
cross-section 
(core) 


Neutron absorption 
cross-section 
(reflector) 


Number of neutrons 
per fission 


Neutron fission 
cross-section 
(core) 


Neutron fission 
cross-section 


Fission energy 


Modified convection 0.0632 cal/em> 


heat transfer 
coefficient 


Reactivity temper- 
ature coefficient 
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TABLE V 


LIST OF COMPUTER RUNS 


RUN METHOD NODES DISTURBANCE COVERGENCE TIME STEP 
CHANGE 
1 CRANK-NICOLSON 45 CENTRAL Ont Tal 
- i 
3 £25 
4 0.01 12 
5 0.001 
6 0.0001 
7 UNIFORM 0.1 
8 SKEWED 
9 72 CENTRAL pe 170 
10 hk 
11 lee 
12 ees 
13 1.4 
14 5 
15 UNIFORM ie 
16 SKEWED 
17 132 CENTRAL 0.1 1.4 
18 a2 
19 15 
20 UNIFORM deed 
21 SKEWED 
a DVOGER 45 CENTRAL Ora: N.A. 
MTH=0 
23 DVOGER OGgt 
MTH=1 
YMAX=1 
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TABLE V (Continued) 


RUN 


24 
25 
26 


27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 


47 


METHOD NODES DISTURBANCE CONVERGENCE 
O01 
001 
DVOGER 45 CENTRAL 0001 
MTH=1 
YMAX=1 
UNIFORM SAE 
SKEWED 
72 CENTRAL 
UNIFORM 
SKEWED 
LZ CENTRAL 
UNIFORM 
SKEWED 
IMPLICIT GEAR 45 CENTRAL dk 
01 
001 
0001 
UNIFORM oak 
SKEWED 
TZ CENTRAL ae 
UNIFORM 
SKEWED 
£32 CENTRAL sil 
UNIFORM 
SKEWED 
45 CENTRAL A) 


46 


TIME..S TEP 
CHANGE 


TART USD 


NAOT THU 
7 = ne Zz 2 
LANTUS 
MAOTTMY 


lawake 


“AOA TMU 


STAG 


=p ts, 


LAR TUS = 


MAOTLAL.. 


CiWaAd 
iA TVD 
MHOIi AU 

CAawl¥2 


MI TVS 


MaO5 LM 


(wine 
LANTHSS 


ad 


VADAAVKOD | ADMARAU Bera 2300"  gOHTHM 


St 


€? AAO TIOLIIMI | 


» 


séf 


TABLE V (Continued) 


RUN METHOD NODES DISTURBANCE CONVERGENCE TIME STEP 
CHANGE 
48 a0 
49 DVOGER 45 CENTRAL Ose 
MTH=1 14 
YMAX=10 
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TABLE VI 
DATA DECK ARRANGEMENT 


CRANK-NICOLSON 
Title 
NUMEL ,NUPBP ,NUMSNP , NFULEL 
List of outer boundary points 
MTH ,MAXDER,NCOUNT 
ZNU, FISFAC,HBAR, EPSVAL, ERRVAL,AFUEL 
TO,H,TF,HMIN,HMAX 
im meucron velocity 
P= “ditftuston “length 
SGA - absorption cross-section 
SGF - fission cross-section 
ALPHA - reactivity temperature coefficient 
PSIIV - initial disturbance flux 
System nodal points with axial and radial position 
Element nodal connectivity 


Nodal neighbor connectivity 


DVOGER 
Title 
NUMEL , NUPBP ,NUMSNP,NFULEL 
List of outer boundary points 
MTH , MAXDER, NCOUNT 
ZNU, FISFAC ,HBAR, EPSVAL, ERRVAL, AFUEL 
TO,H,TF,HMIN , HMAX 
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TABLE VI (Continued) 


ALPHA 
PSIIV 
System nodal points with radial and axial position 


Element nodal connectivity 


IMPLICIT GEAR 
NCORE 
with uniform initial disturbance 
List of core elements 
matle 
NUMEL , NBP ,NUMSNP , NFULEL,NELROW ,MXEVAL,NCOUNT , NCMPK 
NDE ,NL 
List of outer boundary points 
VELOCT , DSUBF , DSUBC, SGMAAF , SGMAAC ,, SGMAF ,, FISFAC ,HBAR 
ZNU , ALPHA, RMSEPS ,AFUEL, TSTART , TEND, PINITV,SCALE 
HMIN , HMAX 
MTH ,MAXDER,NPROB,NTYPE,JPLOT ,NPLOT , IRUN 
Nodal neighbor connectivity 


System nodal points with radial and axial position 


Element nodal connectivity 
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Figure 4. Nodal Neighbor Connectivity 
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